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1. Introduction 

To calculate the hadronic energy spectrum using lattice QCD the procedure is in principle 
straightforward. First one must identify quantum numbers of the channel S of interest which may 
include the particle's intrinsic angular momentum^ J, parity P, charge conjugation C, isospin /, etc. 
Next one constructs the corresponding operators O"^ and O'^^ which will destroy and create a state 
with this symmetry. If one calculates the two-point correlator^ 

G(0 = (o|r<i)^(0<J>^^(o)|o) , (1.1) 

then the energy state spectra must be extracted from the result. Assuming periodic boundary con- 
ditions and a mesonic operator, the fit function is of the form^ 

G(0 = £z„(^'-^"'+^.-^"(^-')) , (1.2) 

where here T is the temporal extent of the lattice. The problem addressed by this paper is how to 
find the coefficients and energy states 

G = {{Zn.En) : n = 1,... , (1.3) 

which minimize x''{G) jndofiG). Here, due to timestep correlations, we have the correlated 
involving the covariance matrix defined by 

X\G) = Y^{G-,-G{mo-%{G-j-G{tj)) , (1.4) 

Oij = G^G] -G'iG'j . (1.5) 
Critically one notes that the number of degrees of freedom, 

ndof{G) = {tmax - tmin + 1) " '^n^ax , (1-6) 

depends on n,„ax, the number of terms in a given fit. Since the latter is unknown, one has a discon- 
tinuous optimization problem having a solution space spanning multiple dimensions. We propose 
the use of an evolutionary algorithm to solve the problem. A complementary discussion of our 
approach may be found in . 



2. Evolutionary Algorithms 

Evolutionary or genetic algorithms use the concept of natural selection to solve function op- 
timization problems. (See and references therein.) The terminology reflects this. Candidate 

' More accurately, since one works on the lattice one is interested in the corresponding lattice quantum number A 
which labels an irreducible representation of the octahedral group which corresponds to the quantum number J of the 
broken continuous rotational symmetry. For instance, 7 = 0— >A = Ai,y = 1 ^A = T]. Higher J correspond often to 
multiple octahedral irreps. See, for example, ^ ^ for mesons and |^, ^] for baryons and references therein. 

^For the sake of the discussion assume translationally invariant zero momentum operators |^ . 

^The asymmetries which arise in baryon correlation functions require only slight modification of this discussion 
which is restricted for simplicity to meson spectra. 
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solutions such as G{t) are organisms. The internal encoding of the solution is its genotype, here 
G = {{Zn,En) :n = \,... ,nmax}- The target function, in our case f{G) = —x'^{G)/ndof{G), is the 
fitness of the organism."* Each step in the algorithm produces a new generation of individuals. 

There are many ways of implementing an evolutionary algorithm. What we use is representa- 
tive: 

1. Create the first generation Pq with randomly generated individuals.^ 

2. Derive P^+i from as follows: 

(a) Mutate each member of the population P^ with a fixed (small) probability.^ 

(b) Select the fittest Neiite organisms and a further Ndiversity random organisms placing all 
their pairwise offspring into Pt+i.^ 

(c) Add N mutant forccd mutations of random elements in the elite to P^+i to explore the 
solution space around the elite. 

3. Repeat until a suitable termination criterion is reached.^ 

In addition to these generic steps one must specify how mutation and breeding are accom- 
plished within the population. For the case of a fit to a single correlator, mutation of an individual 
may include: 

• Adding or removing a random element {Z„,E„) from the genotype's list. 

• Replacing each {Zn,En) by {Zn + l^n:En + where (AZ^jAE,,) are random Gaussian 
deviates.^ 

• Doing a local (e.g. Levenberg-Marquardt) optimization of the fit.*'' 

Breeding or crossover of distinct parent organisms pari and pari produces two child organisms. 
To produce each child, take corresponding ordered pairs in the parent** and generate independent 

^We introduce a minus sign in /(G) to ensure a iiigher value indicates a fitter organism. 
^Here N = {NgUte + N^iversity)^ + N mutant , witli tlie latter constants defined in the algorithm. 
*We do not mutate the fittest organism to ensure that it will survive to the next generation. 

^Here the offspring of the diagonal "pairs" between identical organisms is considered just a copy of the original 
organisms themselves and we thereby are including the elite in the next generation. 

^For instance, one may require a minimum number of generations be exceeded and that a fixed number of genera- 
tions pass with no improvement in fitness of the best organism. A limit on the maximal number of generations may also 
be imposed. 

^Here the standard deviation of the added noise can be tuned to the fitness of our genotype G by making it propor- 
tional to (1 - e-«X'(G)/nrf„/(G)) fpj. ^Qjjjg g^gjj „ 

'"inclusion of such Newtonian optimizations into evolutionary algorithms is often found to be useful We restrict 
the mutation to a fixed number of steps of the local optimization for the sake of efficiency. Greater efficiency might be 
achieved by noting that the linearity of a known fit function of the form ( |l.2| ) admits local optimization with numerical 
methods which exploit the ability to separate linear and non-linear parameters |^] . However such methods will not 
further aid in the discontinuous general problem we are solving of finding said function since here the parameter space 
is not fixed. 

' ' For parents having different numbers of ordered pairs we copy the extra pairs of the longer parent into the child. 
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uniformly distributed random numbers {x,y) E [—5; 1 + 8]}'^ The child element becomes: 

^r^child ^ ^child ^ ^ ^^^parl ^ ( ^ _ ^)2P«'-2 ^ yEP'"'^ + ( 1 - y)EP'"'^) . (2.1) 

Sample fits to single correlators may be found in . 

3. Fitting Multiple Correlators 

Having shown how to fit a single correlator to extract the spectrum it contains, we now gen- 
eralize this to discuss the more practical problem of fitting several correlators. Fitting multiple 
correlators representing a single channel S is desirable as more data allows the resolution of a 
larger number of states with greater accuracy. One typically creates many operators for the channel 
S through group theory methods^^ or by transforming existing operators in ways which conserve 
their channel properties. An example of the former may be seen in figure |l], while quark and link 
smearing of operators which map O'"^ is an example of the latter. The entire set of operators 




Figure 1: Extended gauge field structures contributing the A^'' shown may be combined via octahedral 
Clebsch-Gordan coefficients to quark operators at the origin to create a greater set of operators which couple 
to desired meson channels. The Ty "vector" structure on the right has three components; only the one 
symmetric about the z axis is shown. See [||| and references therein. 

corresponding to S allows one in principle to evaluate an entire correlator matrix Gij{t) between 
them: 

{Of : /=!,..., w}^G,y(0 , (3.1) 

all or some subset of which is to be evaluated and fit. 

Since the correlator matrix grows as the number of operators squared, and because the off- 
diagonal entries have slightly different functional forms, consider the special case of fitting multiple 
diagonal correlators {G;;,/ = 1, . . ./mat} with an evolutionary algorithm. This will require changes 

'^The use of 8 allows for the possibility of extrapolation in addition to interpolation between the parents' parameters, 
thereby avoiding an unwanted rapid contraction to a central point |^ . 

'■'This may be done from the top down by creating an operator space and using group theory projection to extract 
operators with quantum numbers 5 or from the bottom up by using Clebsch-Gordan coefficients to construct the desired 
operators ^ . 
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to our single-correlator problem. First we modify the genotype with {Z,E) — > (Z,/), where index 
/ points to a state list {E\ ,£'m„oJ common to all the correlators.'^ The full genotype becomes 



Fit Genotype = (Dataset coefficients , Energy state list) 

= ((Dataset 1 coefficients, . . .), Energy state list) 

= ((((zj^/j')),...,(z(;j,/;i )),...),(£!,.. .,£,„_)). 

The fitness function /(G) = —x^{G)/ndof{G) is modified due to having multiple datasets to 



(3.2) 



Z'(G) = IZ'(0, 



i=i 



ndof{G) = rid, 



ata f^max 



(3.3) 
(3.4) 



!=1 



where tifiata is the product of the number of timesteps fit and the number of correlators. The com- 
plexity of the genotype permits enhanced evolutionary operations. As a nested hierarchy of lists, 



(3.2) admits more complicated list-based mutations and breeding. Integer indices may be bred and 
mutated bitwise. Finally a reduction mutation which orders masses and coefficients is useful to 
encourage the algorithm to converge to a single representation of the solution. 

In figure ^ the effectiveness of the algorithm to find a known solution is shown. Four synthetic 
correlators, each with 48 timesteps, were created by adding noise to the model function depicted 
on the right of the plot. There one sees four masses displaced horizontally with the corresponding 
coefficients in each respective dataset plotted vertically above the corresponding mass. The left 
side of the plot shows the best fit of each generation. The plot also displays l^dof of the best 
fit (circles) and one sees it converge to 1 as expected as the fit improves.'^ A simultaneous fit to 
actual data of eight diagonal p meson (i.e. A^*' = ) correlators'^ is shown in figure ^. Only the 
energy states are shown of the best fit of each generation up to generation 600. The coefficients in 
each dataset, a further 30 parameters in the final fit, are not shown. The last column depicts the best 
fit found with bootstrap errors produced via Levenberg-Marquardt fits to bootstrap configurations 
with its fixed functional form.'^ 



4. Advantages of Evolutionary Algorithm Fitting 

To conclude, we present advantages of the evolutionary algorithm fitting method. For one it is 
a global optimization method which is furthermore independent of initial conditions. The solution 
space is discontinuous as it spans multiple dimensions, arising from the fact that one does not 
know the exact functional form a priori. The evolutionary algorithm fitting method is capable of 
handling this problem; by minimizing /ridof it finds the number of states in the data in a natural 
manner. As well, the ability to identify whether a state exists or not in an individual correlator in a 

'^Here the integer index / is taken modulo nimax to ensure the coefficient points to an actual energy state. 
'^See [^] for further discussion of this plot. 

'^Simulation details: Wilson quarks, /3 = 6.0, K = .1554, 20"' x 48, 600 configurations, quenched. 

'^Note that as well a Levenberg-Marquardt optimization was done on the best fit found to produce the final result. 
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Figure 2: Masses and coefficients of simultaneous fit to four synthetic correlators. 
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Figure 3: Masses of simultaneous fit to eight p meson (Tj ) correlators. 
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discrete way means that evolutionary fitting can in principle identify wfiich ortliogonal irreps A a 
state straddles and hence aids in the identification of its continuous angular momentum / . 

Computationally, evolutionary algorithms are inherently parallelizable. One can break popu- 
lations into islands breeding on different nodes/CPUs largely independently with only occasional 
migration between them. Large datasets can also be partitioned with sub-genotypes being initially 
evaluated and then stitched together for further evaluation on the entire dataset. Evolutionary fitting 
does not require evaluation of the full correlator matrix, which allows for inexpensive asymmet- 
rical smearing between the sink and source operators.'^ This also means one can restrict oneself 
to evaluating only the diagonal correlators of the correlator matrix where one expects to have the 
strongest signals. This in turn allows a wider assortment of operators to be evaluated. Finally, 
there is the potential for combination with established methods.'^ We are currently applying this 
evolutionary fitting method to analyze all the operators detailed in with promising results. 
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Indeed, our algorithm permits consideration of additional correlators which would not even be allowed in a corre- 
lator matrix, namely those for which only a single operator at the sink or source has the definite quantum numbers of 
channel S. The other operator could potentially couple to many channels beyond that of interest. 

'^One can determine priors for Bayesian fits using distributions of parameters spread across population islands [^. In 
the variational method, one can diagonalize the operators in the usual way and then fit only the new diagonal correlators 
with our fitting algorithm. This would maximize the information in a minimized amount of data for the subsequent 
evolutionary fit. See |p^, for example, on the variational method used with meson correlators. 
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